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Abstract 


A new, nonoscillatory upwind scheme is devel- 
oped for the multidimensional convection equation. 
The scheme consists of an upwind, nonoscillatory in- 
terpolation of data to the surfaces of an intermedi- 
ate finite volume; a characteristic convection of sur- 
face data to a midpoint time level; and a conser- 
vative time integration based on the midpoint rule. 
This procedure results in a convection scheme capa- 
ble of resolving discontinuities neither aligned with, 
nor convected along, grid lines. 

1. Introduction 

The aerospace industry’s wide acceptance of 
computational fluid dynamics is due mainly to the 
successful calculation of transonic flows . 1,3,3 Shock 
capturing schemes, both upwind and central differ- 
enced, have revolutionised the way discontinuous 
flows are investigated. 

In general it is true that an optimally tuned 
and gridded central differenced scheme can produce 
steady state solutions comparable to those produced 
by upwind schemes. It is equally true that upwind 
methods are superior in their ability to resolve un- 
steady flows because central differenced schemes re- 
quire an artificial dissipation term that is usually 
tuned for steady state performance. 

A number of impressive upwind methods have 
been developed for the one-dimensional convection 
equation; many of these schemes have found their 
way into the numerical solution of the Euler and 
Navier-Stokes equations . 4 ’ 6 Most of these schemes, 
with the exception of a few, are second order accu- 
rate in smooth regions and first order accurate at 
extremas. Because this first order behavior can ex- 
cessively damp unsteady calculations, methods that 
are uniformly second order accurate are needed . 3,7 

In general, one-dimensional upwind schemes are 
only formally extended to multidimensions. The 
most common means of extension are the one-step 
Lax-Wendroff, the fractional step, and the multistep 
Runga-Kutta schemes, each with its own benefits 
and drawbacks. The one-step Lax-Wendroff scheme, 
while remaining conservative, requires the evalua- 
tion of cross-derivative terms whose effects on 


accuracy and shock capturing have not been fully 
understood or exploited. The fractional step 
method, by far the most popular means of multi- 
dimensional extension, is not strictly time conserva- 
tive, since all fluxes are not evaluated at the same 
time level. The multistep Runga-Kutta scheme is 
probably the most efficient means of extension but 
unfortunately suffers from an inherent dispersion er- 
ror that can introduce asymmetric behavior. 

The development of true multidimensional up- 
wind schemes has only recently received the at- 
tention previously afforded to formal extensions. 
Davis 8 has developed a rotated scheme that upwinds 
normal to shocks rather than along grid lines, while 
Powell and van Leer 3 have recently formulated a 
convection scheme that obtains its multidimensional 
treatment through a residual distribution step. 

The present approach, in the spirit of van Leer’s 
MUSCL scheme , 13 is to insure nonoscillatory be- 
havior through an uniformly second order accurate 
nonoscillatory interpolation of data to an intermedi- 
ate finite volume. A multidimensional treatment is 
achieved through a characteristic convection of the 
surface data to a midpoint time level. Strict con- 
servation is assured through a midpoint rule tune 
integration. 

2. Analysis 

To illustrate the potential difficulties associ- 
ated with formally extending one-dimensional up- 
wind concepts to higher dimensions, it is sufficient 
to investigate only the one-dimensional linear con- 
vection equation: 


The convection of various gradients (ellipse, 
top hat, and triangle) are solved numerically 
on a uniform grid, using a centered Runga- 
Kutta scheme with artificial dissipation, a one-step 
TVD Lax-Wendroff, and a multistep TVD Lax- 
Wendroff/ Runga-Kutta scheme. Each of the gra- 
dients are constructed over a width of twenty mesh 
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cells and have a nondime nsional amplitude of one 
unit and a length of one half. 

Given the m-stage Runga-Kutta scheme 

u (1) = u n - ai AtR n 
u< 3 > = u" - a 3 Atfl (1) 


i (2) 

« (m > = u" - 

u n+1 = u< m > 

where ax, c*2, etc. are scalar constants and Af is the 
time step. The spatial difference term 

<9u 

J? = — « (u< +1 /2 - Ui_i/ a )/As (3) 

can be constructed from a central difference approxi- 
mation with an added nonlinear artificial dissipation 
term 2 or a flux limited approximation 

«»+ 1/2 = Uj + 0.5&A + u< (4) 

where the convection speed a > 0 is constant; A + 
is a first order forward difference operator; and 4>i is 
the TVD limiter chosen, most often, to take one of 
the following forms 
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corresponding to van Leer’s 11 continous; Roe’s 13 
Minmod and Superbee; and Liou’s 13 exponential 
limiters, respectively. 

These limiters are also used in the construction 
of Sweby’s 14 second-order TVD Lax-Wendroff for- 
mulation 


«? +1 = u" — i/A~Ui — 0.5A <£(l — (6) 

where i/ = a At /Ax and A - is a first order backward 
difference operator. 

To investigate the performance of these 
schemes, a number of test calculations were per- 
formed for eighty time steps at a Courant number of 
one half. The results from a second order, centered, 
four stage Runga-Kutta scheme are clearly oscilla- 
tory (Fig. 1) and may be deemed inappropriate for 
unsteady calculations. It may be conjectured that 


a time dependent artificial dissipation term could 
improve this behavior, but at present no such treat- 
ment exists. 

Results from the one-step TVD Lax-Wendroff 
method (Figs. 2-4) are superior to those obtained 
from the centered scheme and illustrate the relative 
advantages of the upwind method. 

A popular means of extending these TVD con- 
cepts to higher dimensions is through a multistep 
scheme, which unfortunately can introduce undesir- 
able errors. Figures 5 and 6 show that asymmetry 
and oscillation can be introduced even into the one 
dimensional problem and thus a formal extension to 
higher dimensions may be ill advised. 

The asymmetry produced by the Runga-Kutta 
scheme occurs because of the effective limiter cre- 
ated by the multi-step formulation. The use of an 
identical symmetric limiter in each of the multiple 
steps does not insure an overall symmetric behavior. 
The multi-stage formulation has the effect of gener- 
ating nonlinear, spatially shifted terms of the form 
1, which can no longer guarantee a symmetric 
behavior. The use of different limiters in each step 
would be an obvious way of addressing this asymme- 
try and would introduce, in effect, MacCormack’s 16 
predictor-corrector philosophy. 

The multi-step formulation can also degrade, 
with respect to the results from the one-step scheme, 
the performance of the flux limiters. Liou’s expo- 
nential limiter, which performs well in the one-step 
scheme, can be rendered non-TVD and oscillatory 
in the multi-step formulation. This occurs because 
the exponential limiter is TVD for 

r “ (2("“ 1 )/ ri - l) n ^ 

which is unconditionally true only for n = 1. As n 
increases, the limiter has the desirable behavior of 
becoming less diffusive but also produces a more re- 
strictive TVD condition. While this condition seems 
not to be violated in the one-step scheme, it is in 
fact compromised in the first step of the multi-step 
formulation, thus allowing oscillations to develop. It 
again becomes clear that formal extensions to higher 
dimensions may be ill advised. 

Despite the successes achieved in one- 
dimension, the development of the Ultimate’ con- 
vection scheme has proven elusive. Goodman and 
LeVeque 16 have shown that the desire for accuracy 
higher than first order requires sacrificing the TVD 
property in two-dimensions. The obvious next step 
is to retain a non oscillatory behavior despite losing 
the TVD property. This is the rationale taken for 
uniformly second order accurate nonoscillatory one- 
dimensional methods. 

To illustrate the difficulties associated with the 
development of multidimensional schemes, the fol- 
lowing categories are defined: 
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1) Omnidirectionally Nonoscillatory: No new 
extrema* are created in any arbitrary direction. 

2) Preferred Direction Nonoscillatory: No new 
extrema* are created in a preferred direction. 

3) Grid Aligned Nonoscillatory: No new ex- 
trema* are created along at least one grid direction. 


du 


di + a 


du du 

di + b !h 


= o 


a — a(y) and b = b(x) 


( 8 ) 


Most present day nonoscillatory schemes meet the 
requirements of grid aligned nonoscillatory behavior, 
while Davis’s 8 rotated scheme satisfies the definition 
of preferred direction nonoscillation. Definitions 2 
and 3 are necessary but not sufficient conditions for 
an omnidirectional nonoscillatory behavior, which is 
in fact unlikely to be satisfied by anything less than 
the convection of the exact solution. This seem- 
ingly bold statement can be illustrated by a one- 
dimensional convection of a two-dimensional distri- 
bution where dissipation along the convection direc- 
tion can create new extrema* in the cross stream 
direction. A variation of Zalesak’s 17 notched cylin- 
der problem can be used to illustrate this point. A 
notch, fourteen mesh cells wide and twenty mesh 
cell* deep, is cut out of a thirty mesh cell diameter 
cylinder 0.75 unit wide and 0.5 unit high. This dis- 
tribution, Fig.7, can be convected (forty time steps 
at a Courant number of one half, on a uniform grid) 
along the notch’s length-wise direction by a simple 
first-order upwind scheme. The two-dimensional dis- 
tribution uncouple* into a series of one-dimensional 
problems, each of which undergoes some amount of 
amplitude damping. Since the distribution is two- 
dimensional and no information is passed in the cross 
convection direction, the amplitude damping cannot 
be applied uniformly and thus oscillations are cre- 
ated in the cross-convection direction, Fig. 8. As a 
result, it would seem that even a first order TVD 
scheme could be oscillatory in multidimensions, in- 
dicating that the strict nonoscillatory property is 
more restrictive than the TVD property, which is 
the opposite of what is known to be true in one- 
dimension. Based on this observation, one might 
come to the conclusion that, since anything less than 
an exact solution cannot guarantee nonoscillation, 
there is little to be gained by trying to develop a 
truly multidimensional formulation. However, this 
conclusion would be short sighted, since one could 
still make significant improvements in the accuracy 
of unsteady calculations without having to satisfy 
the strict nonoscillatory properties of definition 1. 
In other words, the nonoscillatory property, like the 
TVD property, is not the only measure of accuracy 
and may have to be sacrificed, at least in practice, 
by some other, yet to be established, definition of 
‘goodness’. 

3. Two-Dimensional Formulation 

Consider the two-dimensional linear convection 
equation: 


This equation can be solved numerically, on a 
uniform grid, by a scheme that interpolates grid 
point data to the surfaces of an intermediate finite 
volume; characteristically convects the surface data 
to a midpoint time level; and updates the solution 
with a midpoint rule time integration. 

To insure the production of a nonoscillatory so- 
lution, at least directionally, face values must be con- 
structed in a nonoscillatory manner. To avoid exces- 
sive amounts of damping, this interpolation must be 
at least uniformly second order accurate. 

A second order accurate, nonoscillatory interpo- 
lation, which relies heavily upon a geometric inter- 
pretation of the standard TVD concepts, has been 
recently developed. 7 A one-dimensional interpola- 
tion for a positive convection speed can be written: 

«<+i/2 = + SiAx/2 (9) 

where Si is the slope of a piecewise linear distri- 
bution of data over the intermediate finite volume 
around each grid point. The slope associated with 
the Minmod 1 * scheme can be defined as the median 
of the slopes corresponding to first order upwind- 
ing, central differencing, and second order Warming- 
Beam upwinding. Thus the Minmod slope is 


g Minmod 


where 


A/eA'an(0,u^ 1/2 - u,-,m - u-. 1/2 ) 
Ax/2 


( 10 ) 


«*+l/2 = (11) 

while the slope corresponding to Roe’s Superbee 12 
scheme can be interpreted as: 


gSuperb^e 


0> u ?+l/2 

Median | Ui — u?_ 1//2 , 

Xi j - Uj- 1, Uj+j ~ Uj 

Ax/2 


( 12 ) 


These schemes are second order accurate in smooth 
regions and first order accurate at extrema. To re- 
move this first order behavior, Harten and Osher 6 
developed the uniformly second order accurate 
IJN02 scheme. This scheme can be written as the 
Minmod slope described above, but with a u? +1 ^ 2 
value obtained from a nonoscillatory quadratic in- 
terpolation: 
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u <+i/2 = 0 . 5 (u* + u, + i) - 0 . 25 A+ 1/2 ( 13 ) 

where 

A+i/2 = minmod(Di, A+i) (14) 

Di = D,+i — 2-D* + Di-\ 

and 


of data over a quarter of the intermediate finite vol- 
ume. 

At time level n + 1/2, the surface data is eval- 
uated by first following the characteristic back to 
its spatial location at time level n and then by con- 
structing its value from a bilinear distribution over 
the corresponding quadrant's corner values. These 
quadrant values correspond to the intermediate cell's 
centered, face, and corner values. For positive con- 
vection speeds a and 6, the surface data at time level 
n + 1/2 can be evaluated as: 


minmod(a> b) = 

si0n(a)max(O, «gn(a6)min(|a|, |6|)) ( 15 ) 

This geometric redefinition is rich in extensions 
and has led to the development of Huynh's SONIC 
schemes. 7 By incorporating this new definition of 
u i+i/2 ^ke Superbee slope, one can produce 
Huynh's uniformly second order accurate SONIC- 
B scheme, which is less dissipative than the UNO 2 
scheme* A further extension can be achieved by first 
defining the slope: 


n + l /2 _ n 

“•+1/2,/ ” u <+ 1/2,: 

6 At, 


aAt 




- ^-K> 1/2 .y < + i,x,-x/ 2 ) (2o) 

aJAt 2 , n 

- u"y + u"y_ 1 ( r 3 ) 


where the comer value 


s f = U * +l ^ U *~ 1/2 ( 16 ) 

• Ax v 1 

and then using this value to construct a new slope: 

s sonic-a = Mtdian{0,S\S SONlc - B ) (17) 

This modification results in the SONIC-A scheme, 
which is also uniformly second order accurate but 
less overcompressive than the SO NIC- B scheme. 7 

The median of three quantities, required in the 
sonic interpolations, is evaluated as follows 

medtan(x i, x 3 , X3) = xi 
- hmmmod ( z 3 — zi , Z3 — zi ) ( 18 ) 

while the median of five quantities is evaluated as 


«? +1 /2./-i/a — 0.5(/,/ y -+■ 7y I x ) Ui y (21) 

is evaluated as a symmetric product of I x and J y , 
the one dimensional sonic interpolation operators in 
the x and y directions, respectively. 

Once the surface data has been convected to the 
time level n + 1/2, the finite volume fluxes are eval- 
uated and the values at the grid points are updated 
by the midpoint rule; 


<r = <y - At 


( 


pre + X /3 r . M - 1/2 
■ f »+ l / 2 .j ■*»- 1 /Zj 

Ax 


+ 


^. n + 1/2 ,-, r »+ l /2 

°». j - H /2 ~ °., j - l /2 
Ay 


where F — au and G = 6u. 


( 22 ) 


mcdianixi, x 3 , Z3, Z4, Z5) = medianlXi, X 3 , Z5) 

(19) 

where 

X\ — median(ii, x 3 , x 3 ) 

X 2 — x\ 4 - x 3 -h x 3 + Z 4 

- max(x 1} x 3 ,x 3 ,z 4 ) 

- mtn(zi,X3,z 3 ,z 4 ) - Xx 

The SONIC interpolations are used in each di- 
rection to construct surface data from the grid point 
values at time level n. The surface data is then ad- 
vanced to a midpoint time level using a characteristic 
convection, which is based on a bilinear distribution 


4 . Results 

A series of two-dimensional convection prob- 
lems are used to illustrate the scheme’s ability to 
resolve discontinuities that are neither aligned with, 
nor convected along, grid lines. 

Numerical results are obtained for the convec- 
tion of a cone, box, and cylinder along arbitrary 
flow directions of 0, 45, and 70 degrees. Calcula- 
tions were performed on a uniform grid of 100x100 
mesh cells for 100 time steps at a Courant number 
of one half. Both the cone and cylinder have a 10 
mesh cell radius of 0.25 unit long and are 4 units 
high. The box has a 20 mesh cell square base of 0.5 
unit wide and is also 4 units high. The numerical 
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results using the SONIC-A interpolation, Figures 9 
to 20, show that the various corner, edge, and peak 
discontinuities can be resolved for various flow an- 
gles. The sharp edges associated with the box and 
cylinder distributions are smoothed over a couple of 
mesh cells while the magnitude of the cone’s peak is 
reduced by about six to ten percent. This damping, 
while not unexpected, has been kept to a minimum 
on this relatively coarse grid. These calculations are 
also nonoscillatory in that no trailing, leading, or 
image oscillations are observed in any of these test 
cases. 

Results for the more difficult case of a rotat- 
ing cone and cutout cylinder, 18 Fig. 21, are also 
included to demonstrate the scheme’s shock captur- 
ing capabilities. Calculations were performed on a 
uniform 100x100 mesh cell grid using the Minmod, 
Superbee, SONIC-B, and SONIC-A interpolations. 
Both the cone and the cylinder distributions have a 
30 mesh cell diameter of 0.75 unit wide and an am- 
plitude of 0.5 unit high. The cylinder’s notch is 6 
mesh cells wide and 24 mesh cells deep. The cen- 
ter of rotation is located at the grid center point, 
P(50,50), while the cone is located at P(25,50) and 
the cylinder is centered at P(75,50). These calcula- 
tions were run for one and six complete rotations, 
which correspond to 628 and 3768 time steps, re- 
spectively. This solid body rotation, characterised 
by the velocity field 

a(y) = ~{y - yo) w 

6(x) = (x - xo)w 
x 0 = x(50, 50) (23) 

yo = y(50, 50) 
w = 0.1 

is an extremely difficult test case since the Courant 
number, dispersion, and damping errors vary in 
magnitude throughout the domain, making the sym- 
metric resolution of a symmetric problem highly un- 
likely. 

After one complete rotation, the discontinu- 
ities are captured extremely well by the Superbee, 
SONIC-B, and SONIC-A interpolations. The Min- 
mod results, Fig.22, are extremely dissipative and 
retain only the gross features of the initial distri- 
butions. The Superbee calculation, Fig. 23, resolves 
the notched cylinder well but shows a slight over- 
compression of the cone distribution. The SONIC-B 
interpolation, Fig.24, results in a non-TVD steep- 
ening of the notched cylinder’s narrow bridge but 
has not overcompressed the cone. The SONIC-A in- 
terpolation produces the best overall results, Fig.25, 
since the discontinuities are captured without an ex- 
cessive amount of overcompression or damping. 

The results after six complete revolutions are 
somewhat less encouraging. The Minmod calcula- 
tion, Fig. 26, has all but caused the distributions 


to disappear. The Superbee interpolation, Fig. 27, 
again captures the notched cylinder extremely well 
but compresses the cone into a cylindrical shape. 
The SONIC-B results, Fig. 28, are slightly less over- 
compressive than the Superbee calculation, but are 
asymmetric in form. The SONIC-A results are 
greatly damped and no longer retain the cylinder’s 
plateau region (Fig. 29). However, the cone distribu- 
tion, albeit dissipated, still retains much of its orig- 
inal shape. 

The performance of the SONIC-A interpolation 
is the most desirable from the standpoint of being 
able to retain the overall shapes of the original dis- 
tributions. Since the SONIC-A interpolation is not 
overcompressive, it does not transform gentle gra- 
dients into steep ones, and thus retains the most 
important features of the original distributions. 

5* Concluding Remarks 

A new upwind scheme, based on a uniformly 
second order accurate nonoscillatory interpolation, 
is developed for the multidimensional convection 
equation. Numerical results illustrate the scheme’s 
ability to resolve, without excessive amounts of dis- 
persion or damping errors, discontinuities neither 
aligned with, nor convected along, grid lines. 

It is suggested that the strict nonoscillatory 
property in multidimensions is unlikely to be sat- 
isfied by anything less than the convection of the 
exact solution and is more restrictive than the TVD 
property. If one accepts this belief, it then becomes 
obvious that a more practical definition of ‘goodness’ 
is needed for multidimensional problems. 
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FIGURE 1. - CENTERED 4-STAGE RUNGA-KUTTA CALCULATION. 
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FIGURE 2. - 1-STEP LAX-WENDROFF WITH SUPERBEE LIMITER 
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FIGURE 6. - 2-STEP LAX-WENDROFF WITH LIOU'S EXPONEN 
TIAL LIMITER. N = 10. 



FIGURE 7. - INITIAL CONDITION - 1-D CONVECTION. 
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FIGURE 8. - 1-D CONVECTION OF THE 2-D DISTRIBUTION. 
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FIGURE 9. - INITIAL CONDITION: CONE. 
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FIGURE 10. - 0 DEGREE CONVECTION - CONE. 
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FIGURE 12. - 70 DEGREE CONVECTION - CONE. 
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FIGURE 13. - INITIAL CONDITION - CYLINDER. 




FIGURE 14. - 0 DEGREE CONVECTION - CYLINDER. 
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FIGURE 15. - 45 DEGREE CONVECTION - CYLINDER. 



FIGURE 16. - 70 DEGREE CONVECTION - CYLINDER. 
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MINIMUM = 0.1000 
MAXIMUM = 1.1000 


FIGURE 17. - INITIAL CONDITION - BOX. 
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FIGURE 19. - 43 DEGREE CONVECTION - BOX. 



FIGURE 20. - 70 DEGREE CONVECTION - BOX. 
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MINIMUM = 0.1000 MINIMUM = 0.1000 
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FIGURE 21. - INITIAL CONDITION - SOLID BODY ROTATION. 
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FIGURE 22. - 1 ROTATION - HINHOD INTERPOLATION. 
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FIGURE 23. - 1 ROTATION - SUPERBEE INTERPOLATION. 


MINIMUM - 0.1000 

MAXIMUM = 0.4974 


MINIMUM = 0.0890 

MAXIMUM = 0.5675 



FIGURE 24.-1 ROTATION - SONIC B INTERPOLATION. 
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FIGURE 25. - 1 ROTATION - SONIC A INTERPOLATION. 
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FIGURE 26. - 6 ROTATIONS - HINMOD INTERPOLATION. 
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FIGURE 27.-6 ROTATIONS - SUPERBEE INTERPOLATION. 
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FIGURE 28.-6 ROTATIONS - SONIC B INTERPOLATION. 
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FIGURE 29. - 6 ROTATIONS - SONIC A INTERPOLATION. 


18 


ORIGINAL PAGE 

BLACK AND WHITE PHOTOGRAPH 


Report Documentation Page 

Space Administration 

1. Report No. 

NASA TM- 102354 

2. Government Accession No. 

3. Recipient’s Catalog No. 

4. Title and Subtitle 

A Nonoscillatory, Characteristically Convected, Finite Volume Scheme 
for Multidimensional Convection Problems 

5. Report Date 

6. Performing Organization Code 

7. Author(s) 

Jeffrey W. Yokota and Hung T. Huynh 

8. Performing Organization Report No. 

E-5069 

10. Work Unit No. 
505-62-21 

9. Performing Organization Name and Address 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135-3191 

11. Contract or Grant No. 

13. Type of Report and Period Covered 
Technical Memorandum 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546-0001 

14. Sponsoring Agency Code 

15. Supplementary Notes 

Prepared for the 28th Aerospace Sciences Meeting sponsored by the American Institute of Aeronautics and 
Astronautics, Reno, Nevada, January 8-11, 1990. Jeffrey W. Yokota, Sverdrup Technology, Inc., NASA 
Lewis Research Center Group, Cleveland, Ohio 44135; Hung T. Huynh, NASA Lewis Research Center. 

16. Abstract 

A new, nonoscillatory upwind scheme is developed for the multidimensional convection equation. The scheme 
consists of an upwind, nonoscillatory interpolation of data to the surfaces of an intermediate finite volume; a 
characteristic convection of surface data to a midpoint time level; and a conservative time integration based on 
the midpoint rule. This procedure results in a convection scheme capable of resolving discontinuities neither 
aligned with, nor convected along, grid lines. 

17. Key Words (Suggested by Author(s)) 

Nonoscillatory 

Upwind 

Multidimensional convection 

18. Distribution Statement 

Unclassified — Unlimited 
Subject Category 64 

19. Security Classif. (of this report) 
Unclassified 

20. Security Classif. (of this page) 
Unclassified 

21. No of pages 22. Price* 

18 A03 


nasa form 1626 oct 86 *p or g^g National Technical Information Service, Springfield, Virginia 22161 




